There are potentially many applications that benefit from the use of a computational space that models the chronostratigraphic environment at the time of deposition of the terrains (which is referred to as the “depositional space”). In the following, as a review of the prior art, the specific example of the modeling of the physical properties of geological terrains is put forward in order to provide a clear understanding of the background of the invention.
The physical properties of a 3D geological domain such as porosity or permeability are usually modeled using geostatistical methods. These methods interpolate the ground properties inside a 3D high-resolution mesh on the basis of a usually sparse data, such as measurements obtained from wells. This interpolation process makes an intensive use of distances, such as Euclidean distances between the centers of mesh elements, or of correlation distances given by variograms.
The computed distributions of physical properties must reflect the paleo-environment at the time of deposition of the terrains, so the interpolation between the data obtained from the wells is accurate only if the computed Euclidean distances are close to the equivalent distances at the time of deposition (which are referred to as the “geodesic distances”).
However, rocks have usually been altered since the time of deposition by erosion or faulting and folding caused by applied tectonic stresses. As a consequence, the present-day geometry of a geological domain is usually significantly different from its geometry at the time of deposition.
Thus, geostatistical methods produce inaccurate results if applied directly to the present-day space, described by the Cartesian coordinate system (x, y, z). This problem can be overcome by rather applying geostatistical methods in a “computational space” or “depositional space”, which aims at matching the chronostratigraphic environment at the time of deposition. This computational or depositional space is usually fitted with a curvilinear coordinate system (u, v, w) also called “computational coordinate system” or “depositional coordinate system”.
In all the following, an “isochron” refers to a surface joining points of the present-day space where sediments deposited at the same time. Seismic horizons or the top of geological layers that do not represent a gap in the geologic record are typically isochrons. Boundaries of sequences corresponding to stratigraphic unconformities are not isochrons. “Gaps” refer to situations where points of the depositional space have no corresponding location in the present-day space, and “overlaps” refer to situations where points of the depositional space have more than one corresponding location in the present-day space.
Depositional spaces, as defined in prior art documents, are expected to honor the two following properties:                Property P1: in the depositional space, the geometry of the geological domain represents a chronostratigraphic space where isochrons identified within the stratigraphic sequences are substantially planar and parallel, and where the syn- and post-depositional deformations (i.e. faulting and folding) are substantially removed;        Property P2: every point (u, v, w) of the depositional space located inside a stratigraphic sequence has one and exactly one corresponding location (x, y, z) in the present-day space. In other words, the portion of the depositional space corresponding to a stratigraphic sequence contains neither gaps nor overlaps, allowing points to be mapped from one space to another unambiguously.        
So far, two main types of solutions can be found in the prior art for defining such depositional space.
A first type of solutions can be referred to as “(i, j, k) indexation” solutions, as they consist in building a 3D structured conformal mesh in the present-day space with a Cartesian coordinate system (x, y, z), which represents the present-day geometry of the geological domain. This mesh is structured, so an (i, j, k) index can be assigned to every node and element of the mesh. This indexation is such that the neighbors of a node or an element can be found by simple transformations of its index, these transformations being the same for all the nodes and elements of the mesh. Examples of such transformations are (i−1, j, k), (i, j−1, k), (i, j, k−1), (i+1, j, k), (i, j+1, k) and (i, j, k+1).
Once the (i, j, k) indexation is completed, the depositional space is very simply defined by a curvilinear coordinate system (u, v, w) such as u(x, y, z)=i(x, y, z), v(x, y, z)=j(x, y, z) and w(x, y, z)=k(x, y, z). As the mesh is conformal, there exists a set of facets in the mesh corresponding to every horizon. As a consequence, if all the nodes of these facets have the same k index, the property “P1” is honored. In other words, the mesh must be what is often called a “stratigraphic grid”.
This solution suffers from several drawbacks and limitations:                On one hand, most of the algorithms used for building such stratigraphic grids create badly-shaped (i.e. that is stretched, squeezed, concave or zero-volume) elements in the present-day (x, y, z) space, especially in the neighborhood of faults. As these elements are all right cuboids, or regular parallelepipeds, in the depositional (u, v, w) space, this means that there is a distortion when mapping points from one space to another. In other words, distances measured in the depositional space should be scaled by the distortion to match geodesic distances. This distortion can be very high locally and reduces the accuracy of the property modeling. For instance, it can dramatically alter the computed pore volume of the rocks, and thus the amount of oil and gas in these rocks.        On the other hand, there are some common configurations where it is not possible to ensure that no pair of elements have the same (i, j, k) index, so “null” or “dead” cells are artificially added in stratigraphic grids to solve these problematic configurations. As a consequence, there are gaps or overlaps in the depositional space. Both situations are unacceptable, as they break the required “P2” property.        
A second type of solutions can be referred to as “parametric” solutions, as they consist in computing depositional coordinates with an interpolation of some quantities following a set of geometric rules.
An example of such a geometric solution is given in documents WO 03/050766 by Deny et al. and WO 2005/119304 by Dulac et al., both incorporated herein by reference. Three transfer functions u(x, y, z), v(x, y, z) and w(x, y, z) are interpolated at the nodes of a 3D conformal mesh representing the present-day geometry of the geological domain. The transfer function w(x, y, z) is computed first. In order to meet the requirements described above, the interpolation of the w(x, y, z) transfer function is done so that specific iso-values (iso-surfaces) of the function approximate the geometry of the horizons, which guarantees the property “P1” to be honored. Then, the two other transfer functions u(x, y, z) and v(x, y, z) are computed. The interpolation of these transfer functions is done so that their gradients are orthogonal both to each other and to the gradient of the w(x, y, z) transfer function previously computed, and so that the length of their gradients are equal. Moreover, as an additional boundary condition, the values of these two transfer functions are computed on a reference horizon, specifying that their gradients are orthogonal to each other and that the lengths of their gradients are equal. This reference horizon is chosen so that it intersects a maximum number of faults. WO 2005/119304 discloses running the interpolation on a tetrahedral mesh (a mesh whose elements are all tetrahedra).
This solution suffers from several drawbacks and limitations:                First of all, it falls in the category of geometric and kinematic approaches for modeling the deformation of the terrains through time. As mentioned in the paper “Space-time mathematical framework for sedimentary geology” by Mallet, Mathematical Geology, Vol. 36, No 1, 2004, incorporated herein by reference, computing the curvilinear coordinate system (u, v, w) using the constraints on the gradients which have been shortly described amounts to assuming that terrains have been deformed according to a “flexural slip” kinematic style, that is in which deformation is accommodated by slip along an infinite number of bedding interfaces. As a consequence, the produced depositional space is inappropriate for modeling the physical properties of the geological domain if terrains are believed to deform according to another kinematic style or any other combination of kinematic styles, such as vertical shear, inclined shear, fault-parallel flow, etc. . . . Moreover, it is uncertain which kinematic style or combination of kinematic styles should be used to model the deformation of the terrains. Indeed, kinematic methods are not based on the fundamental principles of the conservation of mass and momentum, which govern rock deformation. In addition, only strain, which is strongly dependent on the kinematic style or the combination of kinematic styles used, can be calculated with kinematic methods, and although Mallet in “Space-time mathematical framework for sedimentary geology” states it is possible to explain how fractures develop and interact, fracture mechanics tells us that the state of stress is rather required.        Secondly, to ensure that points of the depositional space have only one corresponding location in the present-day space, so as to avoid gaps and overlaps and thus breaking the required “P2” property, specific boundary conditions must be applied to the faults when using such a parametric solution: they consist in imposing some arbitrary slip directions and magnitudes for all the nodes of the facets of the mesh matching the faults. This can be considered as a severe limitation, because imposed slip directions and magnitudes are not a priori known and may not correspond to the realistic fault offsets that actually result from the mechanical deformation of rocks. As it is uncertain which slip vectors should be used as boundary conditions, such arbitrary choices can produce a significant distortion when mapping points from one space to another, and as a result the geometry of the geological domain in the depositional space may not be close to its geometry at the time of deposition. In other words, distances measured in the depositional space should be scaled by the distortion to match the geodesic distances. This distortion can be locally very high and reduces the accuracy of the property modeling. For example, it can dramatically alter the computed pore volume of the rocks, and thus the amount of oil and gas in these rocks.        
The purpose of the invention is to propose a method for physically computing a curvilinear coordinate system (u, v, w) defining a depositional space.